% This makes IRFs based on true TFP innovations


eshock = -.01;

% Bivariate
irf2_z = zeros(Tirf,2,varseps2);    %Each page corresponds to different model.
%
for i= 1:varseps2
    xvec2nd = zeros(1,2*lags);      %treat as demeaned.
    if uselagflows ==1 
        %xvec = [4 lags of U/L and 4 lags of another outcome]
        xvec1st = zeros(1,2*lags);
    else
        %xvec = [4 lags of U/L]
        xvec1st = zeros(1,lags);
    end
    epsvec = zeros(1,lags);
    epsvec(1) = eshock;
    %
    for t = 1:Tirf 
        ut= [epsvec, xvec2nd, 0]*ucoef2_z(:,i);
        yt= [epsvec, xvec2nd, 0]*ycoef2_z(:,i);
        xt= [ut;yt];
        %
        xmat = zeros(lags, 2);
        for j = 1:2
            xmat(2:lags,j) = xvec2nd(1+(j-1)*lags : j*lags-1);  % lags=4 and j=1 -> 1:3; j=2 -> 5:7 
            xmat(1,j) = xt(j);
        end
        xvec2nd = reshape(xmat,2*lags,1)';
        if uselagflows == 1
            xvec1st = xvec2nd;
        else
            xvec1st = xvec2nd(1:lags);
        end
        %
        epsvec_ = zeros(1,lags);
        epsvec_(2:lags) = epsvec(1:lags-1);
        epsvec_(1) = 0;
        epsvec = epsvec_;
        %
        irf2_z(t,1,i) = ut;
        irf2_z(t,2,i) = yt;
    end
end
%
%
% Trivariate
irf3_z = zeros(Tirf,3,varseps3);        %Each page corresponds to different model.
%
for i= 1:varseps3
    xvec2nd = zeros(1,3*lags);          %treat as demeaned.
    if uselagflows ==1 
        % xvec = [lags of U/L and f and one other outcome]
        xvec1st = zeros(1,3*lags); 
    else
        % xvec = [lags of U/L and f]
        xvec1st = zeros(1,2*lags);
    end
    epsvec = zeros(1,lags);
    epsvec(1) = eshock;
    %
    for t = 1:Tirf 
        ut= [epsvec, xvec2nd, 0]*ucoef3_z(:,i);
        ft= [epsvec, xvec2nd, 0]*fcoef3_z(:,i);
        yt= [epsvec, xvec2nd, 0]*ycoef3_z(:,i);
        xt= [ut;ft;yt];
        %
        xmat = zeros(lags, 3);
        for j = 1:3
            xmat(2:lags,j) = xvec2nd(1+(j-1)*lags : j*lags-1);  % lags=4 and j=1 -> 1:3; j=2 -> 5:7 
            xmat(1,j) = xt(j);
        end
        xvec2nd = reshape(xmat,3*lags,1)';
        if uselagflows == 1
            xvec1st = xvec2nd;
        else
            xvec1st = xvec2nd(1:2*lags);
        end
        %
        epsvec_ = zeros(1,lags);
        epsvec_(2:lags) = epsvec(1:lags-1);
        epsvec_(1) = 0;
        epsvec = epsvec_;
        %
        irf3_z(t,1,i) = ut;
        irf3_z(t,2,i) = ft;
        irf3_z(t,3,i) = yt;
    end
end
        








